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Abstract 

We consider optimization problems for complex systems in which 
the cost function has a multivalleyed landscape. We introduce a new 
class of dynamical algorithms which, using a suitable annealing pro- 
cedure coupled with a balanced greedy-reluctant strategy drive the 
systems towards the deepest minimum of the cost function. Results 
are presented for the Sherrington-Kirkpatrick model of spin-glasses. 
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1 Introduction. 



There is a standard barrier in applied science: the computational complexity 
of hard (non-polynomial) problems. The modelling of competing interac- 
tions among the components of a large system often lead to consider the 
solution of a problem as the minimum of a functional with a complex land- 
scape. The extensive search for the optimal configurations has a cost that 
grows too quickly (usually exponentially) and become practically intractable 
when the number of composing units is of the order of a few hundreds as 
in the interesting cases. The study of optimizing algorithms is then a basic 
step toward the solution of specific practical problems emerging in different 
fields of applied science. In this paper we build a strategy to efficiently 
explore the landscape of complex functionals in combinatorial optimization 
in order to find its minima both local and global. To allow the reader to 
better focus on our method, let us describe the functional to be minimized 
as the mathematical representation of a quickly changing mountain profile 
(in large dimensions), with a high multiplicity of local minima separated 
by high barriers. The a priori knowledge of the landscape geometry is very 
poor and our strategy to explore the territory in order to find good quality 
minima (close to the global one) is to send signals in random directions (ini- 
tial configurations), follow their evolution according to a specified dynamics 
(algorithm) and collect the observed results. Our investigation procedure 
is not dissimilar from an optical instrument in which we may tune a few 
parameters to better observe the landscape and find the sites which we are 
interested in. The algorithm is preliminary set by choosing the elementary 
dynamical moves: this choice reflects the topology that we are associating 
to our landscape and comes with a notion of vicinity and nearest neigh- 
boring sites. The successive step is to decide the criteria after which to 
select among a large multiplicity of moves. This is done by keeping into 
account what we search for and what we most fear: we want to reach the 
best possible minima as quickly as possible and the worse happening is to 
get stuck in a local minimum which is still far from the optimal or near op- 
timal ones. It appears rather intuitive that an algorithm with a too steepy 
descent (greedy) has a very high risk to get stuck in poor local minima, but 
at the same time a too slow descent (reluctant) would cost a very high price 
in terms of computer time. It is natural to expect, and indeed it is what we 
find, an optimal speed of descent that compromise at best among having 
a wide exploration basin in a reasonable amount of time. Yet the danger 
of remaining caught in wrong local minima remains. To avoid it we also 
allow moves which locally and momentarily deviates from the descending 
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directions. In other terms: to reach a good minimum it is often necessary 
to overcome a high barrier. Physically, the introduction of a similar possi- 
bility works like the availability of thermal energy where the probability of 
its happening is related to the temperature of the system: the higher the 
temperature the more likely are moves upwards and viceversa. To introduce 
such a useful strategy we initially allow upward and downward moves; with 
the time passing the probability to go up is progressively decreased at a 
rate which we may optimize (this simulates the annealing of a physical sys- 
tem) and the algorithm will continue evolving according to its downward 
moves. Our work and the implementation of the algorithm is built and 
tested toward a standard model in combinatorial optimization with origins 
in condensed matter physics: the Sherrington-Kirkpatrick (SK) model for 
the mean field spin glass phase ^1 12] . Among the advantages of our ap- 
proach, there is the flexibility of our algorithms and their wide applicability 
to practical problems like protein folding in biology [2], portfolio optimiza- 
tion in financial mathematics [3], error correcting codes for digital signal 
transmissions |3]. 

2 Results. 

In the following Sections we will present details of the Model and Algorithms 
we used in our simulations. Here we summarize the main ideas and results 
of our analysis. 

In the Sherrington-Kirkpatrick model the cost function is identified with 
the energy of the system, the domain of the cost function is the discrete 
spin configuration space and the optimization problem amounts to find the 
spin configuration with the lowest energy (ground state). Given a proper 
definition of distance in the configuration space (we can think two spin con- 
figurations to be close if they differ only for a single spin- flip), the energy of 
the system is a real-valued function forming a complex and corrugated en- 
ergy landscape, with valleys (local minima) and peaks (local maxima). Our 
optimization algorithms are described as dynamical evolution rules in this 
energy landscape which, starting from a random initial condition, drive the 
system towards local minima of the energy. The random transition from a 
point of the trajectory to the successive, which is a nearest neighboring one, 
is ruled by a probability with exponential density. We consider four differ- 
ent algorithms: starting from the simplest one (Algorithm 0) which allows 
only energy- decreasing trajectories, we implement a sequence of refinements 
(Algorithms 1,2,3) leading to more efficient strategies, which exploit also in- 
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creases in the cost function. 

With Algorithm the cost-decreasing trajectory ends up as soon as it 
reaches a configuration which, according to our notion of vicinity (see Sec. 3) 
is a local minimum. The parameter controlling the transition probability 
function tunes the steepness of descents, generating a continuum of behav- 
iors ranging from a reluctant-type dynamics (very small jumps and slow 
convergence) to a greedy-type one (very large jumps deep into a valley). 

A first improvement of this strategy, implemented in Algorithms 1 and 2, is 
obtained by introducing a "temperature" in the system, which enables ran- 
dom positive fiuctuations of the cost function. This is obtained through the 
choice of a transition probability which gives a non zero weight to upwards 
moves. With these choices we have the following scenario for Algorithms 

1 and 2: the dynamics starts with a given initial temperature and equal 
probability of positive and negative moves. As the time goes on, the system 
is gradually cooled until it reaches a state in which positive fiuctuations are 
forbidden and the dynamics continues as either greedy or reluctant, depend- 
ing on the initial temperature. With a high initial temperature the long term 
behavior of the dynamics will be greedy-like, while a low initial temperature 
will lead to reluctant-type motion. The difference between Algorithm 1 and 

2 lies in the convergence criterium: while the former stops when the first 
local minimum is attained (likewise Algorithm 0), the latter allows the tra- 
jectory to escape from it in view of the possibility to reach deeper minima 
(supplementary stopping conditions are required in this case). 

A further improvement of the algorithm efficiency is obtained with Algo- 
rithm 3. In this case, the transition probability is designed to model an 
initially hot system with high probability of positive moves, which is gradu- 
ally quenched; when the system is cool, positive fiuctuations are absent and 
the decreasing trajectories are forced to follow greedy-like paths. In Fig. ^ 
typical trajectories for the four different algorithms are reported. 

The efficiency of the algorithms are quantified on one hand by measuring 
the average time needed to reach a local minimum, on the other hand by the 
quality of the found minima (i.e. how deep they are). The optimization is 
done by tuning the parameters which control the transition probabilities; in 
particular, for Algorithms 1 and 2 this parameter is mainly the initial tem- 
perature, while for Algorithm 3 it is the rate of the quench, i.e. the speed 
of convergence to zero of the temperature of the system. As one would ex- 
pect, for low initial temperatures (very low possibility of energy increase). 
Algorithm 1 and 2 behaves very much as Algorithm 0. However their differ- 
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ences become effective for sufficiently high initial temperatures. Obviously, 
allowing positive jumps and escapes from local minima, the relaxation times 
increase passing from Algorithm to Algorithm 2; less trivially, numerical 
results show that the scaling of the execution times with respect to the sys- 
tem size is greatly enhanced. This is an important fact, because it suggests 
that a crossover between computation times is to be expected for systems 
with larger sizes. As regards the lowest values found, similar conclusions 
can be drawn: going from Algorithm to Algorithm 2 deeper minima are 
attained. 

Algorithm 3 can be consistently compared with Algorithm 2, which is the 
best performing among the first three. The computation times and their 
scaling with the size are similar for the two algorithms when the initial tem- 
perature (for Algorithm 2) is high, but a clear enhancement is obtained by 
Algorithm 3 when it is low. Also the minimal values of the cost functional 
are similar for high temperatures, while they are lower for Algorithm 2 with 
low initial temperatures. The previous remarks refer to an experimental 
protocol in which the search for low cost configurations is performed testing 
a fixed number of trajectories. The minimization of cost at fixed elapsed 
computer time is another relevant criterium for the comparison of the al- 
gorithms. In this case the best result is obtained with Algorithm 3, even 
though Algorithm 2 gives comparable results. 

3 The model and the algorithms. 

3.1 The Sherrington Kirkpatrick model 

The system we study is the Sherrington-Kirkpatrick model of spin-glasses 
It is defined by the Hamiltonian 



where (Tj = ±1 for i = 1, . . . ,N are Ising spin variables which interact 
through couplings Jij. These are gaussian random variables, independent 
and identically distributed with zero mean and variance 1. The random 
sign (and strength) of the interaction generates frustration in the system, 
i.e. the fact that in low energy configurations some of the couples will have 
unsatisfied interaction. In particular, the ground state of the system is far 
from the standard ground state of ferromagnetic models, where all spins 
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Figure 1: Typical trajectories to reach a local minimum configuration for 
Algorithm 0, 1, 2 and 3. Note that the trajectories generated by Algorithms 
1 and 2 coincide until the first minimum is reached. 

point in the same direction. The model has been solved through the replica 
symmetry breaking ansatz by G. Parisi 0, while the rigorous solution is 
still a debated issue in the mathematical physics community. From the 
numerical point of view, the model poses amazing difficulties and indeed 
it is often presented as the standard example of NP-problems. Several 
numerical studies have tried different algorithms in the search of ground- 
state energies, for example gradient descendent simulated annealing 
13 El, genetic algorithms [HI CDl, extremal optimization [lH [121 ESI- In a 
previous paper we developed a new numerical scheme, which is based on a 
smooth interpolation between greedy and reluctant dynamics [Til El 1^ • 
Here we make a further step by proposing a new class of algorithms which 
we describe in detail in the following. 
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3.2 Dynamical Algorithms 



We focus our attention on stochastic dynamics that generates a sequence 
of spin configurations ending up on a local energy minimum. The smooth 
interpolation between greedy and reluctant dynamics studied in a previous 
work jl6j follows an energy-decreasing trajectory and terminates in the first 
local minimum it encounters: only transitions corresponding to a decrease 
in the cost (energy) function are allowed by the algorithm. In the same 
spirit of Simulated Annealing strategies ^7j, where a slow decrease of the 
temperature leads the system through successive metastable states with 
lower and lower energy, we think of a class of algorithms which also accept, 
in some limited way, transitions corresponding to an increase in the cost 
function. In fact, these algorithms are based on the statistical properties 
of metastable states: they are organized with some structure so that the 
evolution dynamics can be considered as the overlapping of a "fast" motion 
in the basin of attraction of a local minimum and of a "slow" motion with 
jumps between minima (the time of the dynamics is determined by the 
energy barriers between these metastable states). 

In the algorithms that we are going to introduce, the transition between the 
spin configuration at time t, a(t) = (cri(t), . . . , o"Ar(t)), and the successive 
configurations at time t + 1, a{t + 1) = (cri(t + 1), . . . , o"Ar(t + 1)) depends 
on the spectrum of energy changes of (^(t), obtained by fiipping the spin in 
position z, for i = 1, . . . ,N: 



Let also define /S.Ei = mini<j<7v Ai?j that will be used in what follows. As 

a first step, let us briefiy recall the algorithm studied in ^Hl, where only 
energy decreasing trajectory are considered. It is described by the following 
procedure: 

Algorithm 

1. Initialization: choose an initial spin configuration a{0) and a param- 
eter value for A > 0. 

2. Generate a random number D with probability density 




(2) 




(3) 
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3. Select the site i* associated with the closest energy change to the value 
D, i.e.: 

i* : \AEi* -D\= min {lA^;^ - D\ : AEi < 0}. (4) 

ie{i,...,N} 

4. Flip the spin on site i*: 

/, i\ f —(^iit) if i = i* 

+ = if>#>.. (5) 

5. If AEi > 0, Vi = 1, . . . , A/", then the algorithm stops {cr{t) is a local 
minimum); otherwise repeat from step 2. 

The dynamics generated by this algorithm follows a 1-spin flip decreasing 
energy trajectory and arrives at a configuration whose energy cannot be 
decreased by a single spin-flip. The control parameter A in the probability 
distribution ftmction for the move acceptance, tunes the speed of conver- 
gence to local energy minima: the larger is A, the bigger is the probability 
of doing small energy-decreasing steps, so that the trajectory will follow 
an evolution path close to level curves (reluctant) while, small values of A 
enrich the probability of large negative energy steps (greedy), which will 
quickly drive the dynamics to the end-point. 

As a modification of Algorithm we consider two new algorithms (Algo- 
rithm 1 and Algorithm 2). They generate a dynamics that follows a 1-spin 
flip trajectory that, in addition to energy-decreasing transitions, accepts 
also energy- increasing transitions with probability exponentially decreasing 
in time. The difference between the two is that while the trajectory of Al- 
gorithm 1 ends up in the first local minimum it encounters, in Algorithm 2 
it may continue to explore the space of configurations through the visit of 
subsequent local minima. 

Algorithm 1 

1. Initialization: choose an initial spin configuration a(0) and parameter 
values < ci(0) < Ai, < C2 < A2(0), with the obvious constraint 

In our simulation we chose Ai as the only free parameter, by taking 
A2(0) = Ai, ci(0) = Ai/2, C2 = Ai/2. This amounts to start with an 
equal probability of energy decreasing and energy increasing transi- 
tions (ci(0)/Ai = C2/A2(0) = 1/2). 
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2. Generate a random number D with probability function 

ci(t)e^i^' ifx<0 
C2e 



3. Select the site i* associated with the closest energy change to the value 
D and with the same sign, i.e.: 

i* : \^Ei* -D\= min {lAE, - D\ : AE, ■ D > 0}. (8) 

ie{i,...,N} 

4. Flip the spin on site i*: 

cr^it + 1) = I ^J^ = t (9) 

5. If AEi > 0, Vi = 1, . . . , A^, then the algorithm stops {<j(t) is a local 
minimum). Otherwise, change the parameter X2{t) of the probability 
distribution in step 2 with a suitable scheduling, for example 

A2(t) = < A; < 1 (10) 

and return to step 2. 

The trajectory generated by Algorithm 1 wonder in the energy landscape 
(by a succession of moves which decrease and increase energy) till it arrives 
to a local minimum. Starting from a symmetric probability distribution for 
the spin-flip selection, as time goes on the probability of energy-increasing 
moves is decreased by the update rule ([TU|]. 

Next, we want to consider an algorithm as the previous one but with the 
possibility of exploring subsequent minima. The problem one has to solve 
is to give an efficient criterium to stop the dynamics. We considered the 
following implementation: 

Algorithm 2 

1. Initialization: as in Algorithm 1. Set also m = 1000 and e = 10""^. 

2. Generate a random number D as follows: 
with probability function 

/ ci(t)e^i- if x<0 ci(t) C2 

= 1 c,e-^«^ if . > ^ ^ 
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and with probability function 



J[x) = < „ „ if — — > m— — - (12) 

^ ^ \ if X > Ai X2{t) ^ ' 

3. Select the site i* associated with the closest energy change to the value 
D and with the same sign, i.e.: 

i* : \/\Ei*-D\= min {I AE^ - : AE^ ■ > 0}. (13) 

i&{l,...,N} 

4. Flip the spin on site i*: 

5. If AEi > 0, Vi = 1, . . . , N, and Pt(/^ > AEj) < e then Stop. 

D is a random number, Pt is the cumulative function of the probability 
described in step 2 and e is a small parameter. In other words, if 
we arrive in a minimum and the probability of a significant energy 
increasing transition from this local minimum is too small (or even 
zero when the energy increases are forbidden, see step 2), then the 
algorithm stops. 

6. Change the probability distribution (jll|) with the scheduling (jlUj) for 
A2(t) (the same scheduling used in Algorithm 1) and return to step 2. 

As in Algorithm 1, the dynamics generated by this algorithm follows a 1- 
spin flip trajectory making a combination of upwards and downwards moves. 
However, in this case, the trajectory does not end up in the first 1-spin flip 
stable configuration it encounters, at least as long as the probability of 
positive moves (c2/A2(t)) remains greater than a certain threshold (1/m 
times the probability of negative moves ci (t) / Ai - in our experiments m = 
1000). With this strategy it is possible to escape from the local minima to 
explore the neighboring space in view of (possible) lower energy minima. 
When the probability of energy increases exceed this fixed threshold, from 
this point on, only decreases in energy are accepted and so the process 
terminates when the subsequent local minimum is reached. In fact, when 
the process starts at time t = we choose equal probabilities ci(0)/Ai and 
C2/A2(0) of cost-decreasing or cost-increasing moves, respectively, by settling 
C2 = A2(0)/2. As the algorithm continues its execution, we decrease C2/A2(t) 
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towards zero, varying the control parameter X2{t) in accordance with the 
above mentioned law (jlUj) : 

A2(t) = ^, A2(0) = Ai, 0<A;<1 

(and keeping fixed Ai) until < "m-j^', as a consequence, the probability 
of energy- decreasing move acceptance Ci(t)/\i tends to one (ci(t) = Ai(l — 
C2/A2(t))). Therefore, while the speed of convergence to the local energy 
minima is mainly tuned by Ai, the vanishing velocity of the probability of 
energy-increasing steps is governed by the parameter k. Of course, large Ai 
(and A2(t)) lead to evolution paths generated by small (in absolute value) 
energy changes [annealed reluctant dynamics) and the closer A; is to 1, the 
slower X2{t) grows and then the more energy increases are enabled. When 
> ''TT'j^ the dynamics continues governed only by the parameter Ai, 
not depending on t. 

We see that for Algorithm 2 the possibility to escape from the minima is 
effective only when Ai is sufficiently small (say Ai ~ 1, and then A2(0) ~ 1, 
see (Uni))- For greater values of Ai the possibility to explore successive 
minima is not exploited and both the dynamics 1 and 2 can be expected 
to give similar results in terms of achieved minimum energy level. In these 
cases, the dynamics generated by Algorithm 2 ends up naturally, after t' 
steps, in the first minimum it encounters, because the (step dependent) 
probability Pt' to escape from this configuration is too small; therefore, we 
expect that for large values of Ai Algorithms 1 and 2 should be equivalent. 

Since for these algorithms the speed of convergence to the finale state is 
governed by the probability function ft{x), we can consider a third algorithm 
in which the time dependence is present only in the control parameters Aj(t), 
i = 1,2] in this case, starting from a (in general) non symmetric probability 
function, the dynamics evolves gradually towards a final scenario in which 
the system is cooled by tuning the control parameter Ai(t). 

Algorithm 3 

1. Initialization: choose an initial spin configuration a{0) and parameter 
values Ai(0), A2(0) such that 1/Ai(0) + 1/A2(0) = 1. Set also m = 1000 
and e = 10"^ 

2. Generate a random number D as follows: 
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with probability function 



^ e -^^Wx if a; > Ai(t) " X2{t) 

and with probabihty function 

. / Aie^i^ ifx<0 1 ^ 1 

fix) = S T r> II , , , > m , , , 16) 

^ ^ \ if X > Ai(t) A2(t) 

3. Select the site i* associated with the closest energy change to the value 
D and with the same sign, i.e.: 

i* : \AEi*-D\= min {\AEi - D\ : AEi ■ D > 0}. (17) 

ie{i,...,Af} 



4. Flip the spin on site i*: 



—(Ti{t) if z = i* 
crj(t) if i 7^ 



5. If AEi > 0, Vi = 1, . . . , iV, and Pt{D > AEi) < e then Stop (as in 
Algorithm 2) . 

6. Change the probability distribution defined in ()15|) with the same 
scheduling for \2(t) used in Algorithm 2 and return to Step 2. 

The main difference between Algorithm 2 and Algorithm 3 is that in the 
latter, when the process starts at time t = we have (if Ai(0) ^ 2) dif- 
ferent probabilities of energy- decreasing moves (1/Ai(0)) and of energy- 
increasing moves (1/A2(0)). As Algorithm 3 continues its execution, we 
decrease 1/A2(t) towards zero, varying the control parameter A2(t) in accor- 
dance with the scheduling: 

Mi) = ^, VO) = ^Aa_. o<t<l (19) 

until < rrij^; as a consequence, the probability of energy-decreasing 

move acceptance 1/Ai(t) tends to one (Ai(t) = -^^^^zj)- Therefore, while the 
speed of convergence to the final state is mainly tuned by the initial value 
Ai(0) of the time dependent parameter Xi{t) (which tends to 1, as time 
t increases), the vanishing velocity of the probability of energy-increasing 
steps is governed by the parameter k. When xr(t*) ^ ''^x^F) dynamics 
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Figure 2: Probability density functions for Algorithm 1 and 2 (part (a)) and 
for Algorithm 3 (part (b)) for different values of time t. The continuous lines 
refer to t — 0; the time goes on passing from broken hues to dotted ones. 
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continues, for t > t*, governed only by the parameter Ai = Ai(t*) (close to 
1) not depending on t. The dynamic evolution of the probability density 
functions for Algorithm 1 and 2 compared with Algorithm 3 is reported in 
Fig. El 

Summarizing: the control parameters are A for Algorithm 0, Ai and k for 
Algorithms 1 and 2, and Ai(0) and k for Algorithm 3. Varying them we 
study the efficiency of the algorithms by measuring the average time to reach 
a metastable configuration and the lowest energy value found for different 
system sizes. 

4 Data analysis. 

To compare these annealed algorithms with those carried out in previous 
works ^JEl^] and in particular with Algorithm 0, we performed a set of 
trials for different values of A^, starting from initial conditions (for a sys- 
tem of size A^) and averaging the data on ureal = 1000 disorder realizations. 
We measured two quantities to test the performance of the algorithms: 

- the average time (i.e. the number of spin fiips) to reach a minimum 
energy level 

1 ^ 

r=^E^- (20) 

i=l 

with M = N ■ ureal and ti, i = 1, . . . , M the time for each initial 
condition; 

- the lowest energy found (averaged over disorder) 

/ min^ij'jv(J,a) \ 

where min^ H]si{J, a) is the minimum value of the energy of the meta- 
stable states attained starting from the set of the A^ initial conditions. 

Our numerical experiments follows two different protocols: 

1. with a fixed number of initial conditions; 

2. with a fixed elapsed computer time. 

The results are described in the following subsections. 
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Figure 3: Average time r to reach a metastable configuration as a function 
of N for different values of Ai and k for Algorithm 2 and for a fixed number 
of initial spin configurations. 
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4.1 Fixed number of initial conditions 

The dynamics of Algorithm has been shown UB] to behave as a smooth 
interpolation between greedy and reluctant dynamics J3] depending on the 
parameter A: small A (say A ~ 1) plays the role of the greedy algorithm, 
while large A (say A ~ 100) that of reluctant. In fact, the relaxation time 
t{N) grows linearly with the system size when A ~ 1 and quadratically when 
A ~ 100 (see Tab. H}, as it was previously observed in [T3] for deterministic 
greedy and reluctant regimes. 

In Fig. 01 which refers to Algorithm 2, we represent r as a function of 
(A^ G [25, 300]). We performed the analysis for different values of the control 
parameters. For the sake of space, we show only the values Ai = 1, 10, 100 
and three values of k {k = .98, .99, .995) for each Ai, together with the best 
numerical fits. Fig. El shows the progressive increase of the slope in log-log 
scale from a sub-linear law in A^ for Ai = 1 and k = .98 (-^) to a super-linear 
one for Ai = 100 and k = .98 (■ x-). More in detail, the numerical fits of 
T\i,k{N) ~ A^" in Fig. 121 are reported in Tab^ 



Table 1: Numerical fits of tx{N) ~ A^'' for Algorithm (with the symbols 
of Fig. ini) and of TA^,fc(A^) ~ A^"^ for Algorithm 1 and Algorithm 2 (with the 
symbols of Fig. ^ 
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With the same protocol (fixed number of initial conditions), we measured 
the lowest energy Hj^j found by the algorithms. As a general remark we 
recall that from a theoretical point of view it is proved the monotonicity in 
A" of the ground state energy (this follows from sub-additivity [l^)- For the 
largest size we have studied, some values of the simulation parameters give a 
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non-monotone behavior in N, suggesting that we are not actually finding the 
true lowest energy state. A larger number of trials (i.e. initial conditions) 
would be needed to achieve the global minimum. However, our principal 
aim here is not to have a perfect measure of ground state energies. In Fig. |3] 
we represent, for Algorithm 2, as a function of N for different values of 
Ai and k. The best results for large N are obtained for Ai = 100 and k = .98 
which corresponds to annealed reluctant dynamics (as found for Algorithm 
0, see Fig. EJ. Therefore, this confirms [121^21 that, for a fixed number of 
initial spin configurations, the algorithm that makes moves corresponding 
to the "smallest" possible energy change keeping the possibility of energy 
increase only for the first steps of the algorithm is the most efficient in 
reaching low-energy states. Note that, for Ai = 1 and k = .995 the attained 
energy values are sufficiently low: even if these results are not better than 
those for Ai = 100 (with k = .98 and k = .995), they should not be 
discarded since the average time scales better ( Ti^995(A^) ~ A^-^^^ instead 

of rS,.98(iV) - N'-''' or 4t.995(A^) ~ N^'''') 

Comparing these results with those obtained with the interpolating 

greedy and reluctant algorithm (Algorithm 0) we note (Figs. and 

|B]and Tab. that for small A and Ai Algorithm 2 is better performing than 

Algorithm both with respect to average time and energy levels, while 

for greater A and Ai we find comparable energy values but with lower cost 

for the computational time for Algorithm 2 (t^qq 98(A^) ~ A^^ '''^^ instead of 

The same analysis is considered also for Algorithm 1. The comparison 
between Algorithms 1 and 2 shows that the possibility of exceed the energy 
barriers between minima is useful only for small values of Ai (for Ai close 
to 1 Algorithm 2 is more efficient than Algorithm 1 in reaching lower en- 
ergy states) while for Ai > 5 the performances of Algorithms 1 and 2 are 
practically indistinguishable (see Figs. [7| and |H1). Moreover, we note that 
the best scaling of the average time Tx-^^k with respect to is obtained with 
Algorithm 2 (see Tab.[T)), though for fixed A^, Ai and k, we have t^^^ < r^^^,- 

Figures IHl and ^1 report the results of the analysis of Algorithm 3 with a 
fixed number of initial conditions: A^ G [25,400]) for three distinct values 
of Ai(0) (Ai(0) = 2, 10, 100) and for four values of k {k = .98, .99, .995, .997) 
for each Ai(0). Because of high computational costs (which increase with 
Ai(0) and k), the cases A^ = 350 and A^ = 400 for Ai(0) = 100 are only 
partially studied. For the same reason also the case k = .997 is considered 



"'^From now on, the superscript (x) in the notation of the average time r^^-* will refer 
to the number of the corresponding algorithm. 
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gure 4: Lowest energy value H^^ as a function of for different values of 
and k for Algoritlim 2 and for a fixed number of initial conditions. 
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Figure 5: Lowest energy value Hn as a function of obtained using a 
protocol with a fixed number of initial conditions for A = 1 (*) and A = 100 
{<>) for Algorithm and for Ai = 1 and k = .995 (□) and for Ai = 100 and 
A; = .98 (x) for Algorithm 2. 
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Figure 6: Average time r to reach a metastablc configuration as a function 
of for A = 1 (+) and for A = 100 (o) for Algorithm 0, and for Ai = 1 and 
k = .995 (□), for Ai = 100 and k = .98 (x),and for Ai = 100 and k = .995 
( *) for Algorithm 2. 
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Figure 7: Lowest energy value as a function of N for Ai = 1 and different 
values of /c, for Algorithm 1 and 2. 
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Figure 8: Lowest energy value as a function of N for Ai = 10 and for 
different values of k obtained with Algorithm 1 and 2. 
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only for Ai(0) = 2. 

Fig. ^1 shows that Algorithm 3 seems to depend weakly on the parameter 
Ai(0), its behavior being mainly ruled by k. In fact, the lines of the if^v 
values corresponding to the same choices of k are grouped into narrow bands 
well separated one from the others. Moreover, a closer look to Fig. 1101 shows 
that the best result for Hn is obtained for Ai(0) = 2 and k = .997. Note 
that for any Ai(0), the closer the values of k to one, the lower the values 
of energy: slow growths of the parameter X2{t) enable energy increases and 
then the possibility to exceed the energy barriers. Even though Algorithm 
2 is slightly better performing (Ai = 100, k = .98 see Fig. 0} in terms of 
minimum energy level reached, the best scaling of rAi(o),fc(A^) is obtained 
by Algorithm 3. In fact, for Algorithm 3 we note (Fig. El and Tab. 
the progressive increase of the slope in log-log scale from a scaling law 
^xX),kiN) ~ for Ai(0) = 100 and k = .995 (■*■) to r^^^^^AN) ~ N-'^ 
for Ai(0) = 2 and k = .98 ( ^). More in detail, the numerical fits of 
''"Ai(o),A;(-^) ~ for Algorithm 3 are reported in Tab|21 

To conclude the analysis of the protocol with a fixed number of initial 
conditions we can say that taking into account also the average time r, the 
best performing algorithm in reaching minimum energy level is Algorithm 3 
(Fig. Ill|) . In fact. Algorithm 3 with Ai(0) = 2 e k = .997 attains minimum 
energy levels comparable with those obtained by the other algorithms with 
A and Ai equal to 100 but with lower computational costs (rg^ggy ~ N-'^'^'^ 
while rSJ ~ ^« ^ ivi.724 g^^^ ^(2)^ ^^ _ ^1.771^ Tabs. O and 

121. 

4.2 Fixed elapsed computer time 

Finally, we analyze the lowest energy states found by the dynamics varying 
the control parameters for a given elapsed running time for all algorithms. 
In Fig.^Jwe consider the minimum energy values Hn, obtained by choosing 
different system sizes and, for each of them, different parameter values 
(A = 1,10,100 for Algorithm 0, Ai = 1,5,10,100 for Algorithm 1, Ai = 
1, 10 for Algorithm 2 and Ai(0) = 2, 10, 100 for Algorithm 3) with different 
annealing scheduling each {k = .98 and k = .995 for Algorithm 1 and 2, 
k = .995 and k = .997 for Algorithm 3), for a fixed time of 50 h of CPU on a 
IBM SP4. For Algorithm 2 we consider in detail mainly the case (Ai = 1) in 
which the dynamics behaves differently from that generated by Algorithm 1 . 
Each run (i.e. for fixed N and for fixed control parameter) consists of 1000 
disorder realizations, with the same CPU time length (3 min.) assigned to 
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each sample, in order to compare these results with ^1^1^]. With all this 
dynamics, for < 150, we believe to find the ground state of the system, 
since varying the control parameters and independently on the algorithm 
used, the values of coincide, within our numerical accuracy (10~^°). 
The best result is obtained with Algorithm 3 for the case Ai(0) = 2 and 
k = .997 (even though the result provided by Algorithm 2 for Ai = 1 and 
k = .995 is comparable). Note that, for Algorithm 1 the best result is for 
Ai = 10 and /c = .98 in good agreement with the best result of Algorithm 
obtained for A = 10 fFig. I12|). Moreover, it is worthnoting that the values 
Hjsr obtained with Algorithm 3 for the case Ai(0) = 2 and k = .997 are the 
best (for fixed CPU time) with respect to all algorithms we consider in the 
present paper and in jTH [1^1 CHI ■ 



Table 2: Numerical fits of ta^(o),a:(^) ~ N"- for Algorithm 3 
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Figure 9: Average time r to reach a metastable configuration as a function 
of for different values of Ai(0) and k for Algorithm 3, together with the 
best numerical fits for a fixed number of initial conditions. We represent 
Ai(0) = 2{k = .98 (^), k = .99 (+) and k = .995 (-B)), Ai(0) = 10 (A; = .98 
(<3>-), A; = .99 (-f •) and k = .995 (•□•)) and Ai(0) = 100 {k = .98 (-x-), k = .99 
(•A-) and A; = .995 (*•)) 
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Figure 10: Lowest energy value Hj^ as a function of N for different values 
of Ai(0) and k for Algorithm 3 and for a fixed number of initial conditions. 



26 



-0.67 




-0.75 I 1 1 1 1 1 1 1 

50 100 150 200 250 N 300 

Figure 11: Lowest energy value as a function of for A = 100 (o) for 
Algorithm 0, for Ai = 100 and k = .98 (A) for Algorithm 1 and (x) for 
Algorithm 2 and for Ai(0) = 2 and k — .997 (*) for Algorithm 3. 
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Figure 12: Lowest energy value Hjq as a function of N for different values 
of control parameters for Algorithms 0, 1, 2 and 3, for a fixed CPU time of 
50 h on a IBM SP4. The symbol (+) refers to A = 10 for Algorithm 0, (o) 
to Ai = 10 and k = .98 for Algorithm 1, (□) to Ai = 1 and k = .995 for 
Algorithm 2 and (x) for Ai(0) = 2 and k = .997 for Algorithm 3. 
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